%CODE FOR PART2 OF SS 2015;
close all; clear; clc;


tic;
%% Start Here: Setup working folder, seed, and options;
    %n Set Seed;
state=1991;
rand('state',state);
randn('state',state);
stream1= RandStream('mt19937ar','Seed', 1);
RandStream.setGlobalStream(stream1);

%% 0. Import Data and Definitions;
    
    %Import Extended DW Nominate scores;
data_smdwnom=importdata('dwnominate_sm2015.csv',',',1);

dw_congress_ext=cellfun(@str2num,data_smdwnom.textdata(2:length(data_smdwnom.textdata),1));
dw_state_ext=cellfun(@str2num,data_smdwnom.textdata(2:length(data_smdwnom.textdata),3));

dw_extended_112=data_smdwnom.data(:,1);
dw_extended_111=data_smdwnom.data(:,2);
dw_extended_se=data_smdwnom.data(:,5);

    %Import original regression data from S&M (2015)- at state level;
N=110;
data_sm=importdata('data_sm2015.csv');

sm_congress=cellfun(@str2num,data_sm.textdata(2:length(data_sm.textdata),4));
sm_state=cellfun(@str2num,data_sm.textdata(2:length(data_sm.textdata),5));
sm_districts=cellfun(@str2num,data_sm.textdata(2:length(data_sm.textdata),15));

    % Create variables needed for regression;
libs_smorig=cellfun(@str2num,data_sm.textdata(2:length(data_sm.textdata),37));
c=ones(N,1);
stfracmm= cellfun(@str2num,data_sm.textdata(2:length(data_sm.textdata),16))./sm_districts;

demprev=cellfun(@str2num,data_sm.textdata(2:length(data_sm.textdata),6));
repprev=cellfun(@str2num,data_sm.textdata(2:length(data_sm.textdata),7));
demnext=cellfun(@str2num,data_sm.textdata(2:length(data_sm.textdata),8));
repnext=cellfun(@str2num,data_sm.textdata(2:length(data_sm.textdata),9));

natdprev=cellfun(@str2num,data_sm.textdata(2:length(data_sm.textdata),10));
natrprev=cellfun(@str2num,data_sm.textdata(2:length(data_sm.textdata),11));
natdnext=cellfun(@str2num,data_sm.textdata(2:length(data_sm.textdata),12));
natrnext=cellfun(@str2num,data_sm.textdata(2:length(data_sm.textdata),13));

    % Here are the two votlib2 definitions, we use the second one, which is
    % the one used by S&M (2015);
%votlib2=((demprev+demnext)./(demprev+demnext+repprev+repnext))-((natdprev+natdnext)./(natdprev+natdnext+natrprev+natrnext));
votlib2=demnext-natdnext;

repgerr=data_sm.data(:,3);
bipartgerr=data_sm.data(:,4);

blackpop=cellfun(@str2num,data_sm.textdata(2:length(data_sm.textdata),31));
hispop=cellfun(@str2num,data_sm.textdata(2:length(data_sm.textdata),32));

interac_rep=stfracmm.*repgerr;
interac_bi=stfracmm.* bipartgerr;

libslag_smorig=cellfun(@str2num,data_sm.textdata(2:length(data_sm.textdata),38));
  
%% Part 2.1. Replicate SM regression and do Monte Carlos simulation;

    %Replicated Conditional Effects Model Variables;
[libs_sm,libslag_sm]=create_liberal(sm_congress,sm_state,dw_congress_ext,dw_state_ext,dw_extended_111);


Y_sm=libs_sm;
X_sm=[stfracmm, repgerr,bipartgerr, interac_rep,interac_bi, votlib2,blackpop,hispop,libslag_sm,c];
w_sm=sm_districts;

    %Weighted OLS estimates and standard errors from Conditional Effects Model;
[b_sm,se_sm,tstat_sm]=ols_weighted(X_sm,Y_sm,w_sm)


    %Monte Carlo parametric procedure: for each member draw noise from normal with changing levels of standard deviation,
    % procedure as in first part of paper ;
rep=1000;
delta=0.01;
error_max=0.2;
grid=0:delta:error_max;

bsm_mc=cell(1,5);
ssm_mc=cell(1,5);
tstatsm_mc=cell(1,5);

for i=1:5
    bsm_mc{i}=zeros(rep,length(grid));
    ssm_mc{i}=zeros(rep,length(grid));
    tstatsm_mc{i}=zeros(rep,length(grid));
end
    %Original results- noise level zero;
for i=1:5
    bsm_mc{i}(:,1)=b_sm(i).*ones(rep,1);
    ssm_mc{i}(:,1)=se_sm(i).*ones(rep,1);
    tstatsm_mc{i}(:,1)=tstat_sm(i).*ones(rep,1);
end

    
for i=2:length(grid);
    %Here I draw the normal noise given standard deviaiton grid(i);
    noise=normrnd(0,grid(i),length(dw_extended_111),rep);
    for j=1:rep;
        dw_extended_noise=dw_extended_111+noise(:,j);
        
        [libs_noise,libslag_noise]=create_liberal(sm_congress,sm_state,dw_congress_ext,dw_state_ext,dw_extended_noise);
    
        Y_noise=libs_noise;
        X_noise=[stfracmm, repgerr,bipartgerr,interac_rep,interac_bi,votlib2,blackpop,hispop,libslag_noise,c];

        [b_noise,se_noise,tstat_noise]=ols_weighted(X_noise,Y_noise,w_sm);
        
        for l=1:5;
            bsm_mc{l}(j,i)=b_noise(l);
            ssm_mc{l}(j,i)=se_noise(l);
            tstatsm_mc{l}(j,i)=tstat_noise(l);
        end
    end
end

    % Plot parametric results;
t_critical=norminv(0.975,0,1);

b_pcpossig=zeros(length(grid),5);
b_pcnegsig=zeros(length(grid),5);

    % Significaint coefficients;
for i=1:5;
    b_pcpossig(:,i)=mean(tstatsm_mc{i}> t_critical)';
    b_pcnegsig(:,i)=mean(tstatsm_mc{i}< -t_critical)';
end
    % 10 graphs, share positive and negative coefficients for each level of
    % noise, for the five variables of interest;
fontsize=14;
width_line=4;
ylim_up=1;

        % B1-Fraction Majority Minority;
figure('units','normalized','position',[.1 .1 0.5 0.2]);
subplot(1,2,1);
plot(grid, b_pcpossig(:,1),'k', 'LineWidth',width_line);
ylabel('Proportion Positive and Significant', 'FontSize',fontsize)
ylim([0 ylim_up])
xlabel('Standard Deviation of Error', 'FontSize',fontsize)
xlim([0 error_max])


subplot(1,2,2);
plot(grid, b_pcnegsig(:,1),'k', 'LineWidth',width_line);
ylabel('Proportion Negative and Significant', 'FontSize',fontsize)
ylim([0 ylim_up])
xlabel('Standard Deviation of Error', 'FontSize',fontsize)
xlim([0 error_max])

set(gcf,'PaperPositionMode','auto');
saveas(gcf, 'figure2_1a.png');

        % B2-Republican Control;
figure('units','normalized','position',[.1 .1 0.5 0.2]);
subplot(1,2,1);
plot(grid, b_pcpossig(:,2),'k', 'LineWidth',width_line);
ylabel('Proportion Positive and Significant', 'FontSize',fontsize)
ylim([0 ylim_up])
xlabel('Standard Deviation of Error', 'FontSize',fontsize)
xlim([0 error_max])

subplot(1,2,2);
plot(grid, b_pcnegsig(:,2),'k', 'LineWidth',width_line);
ylabel('Proportion Negative and Significant', 'FontSize',fontsize)
ylim([0 ylim_up])
xlabel('Standard Deviation of Error', 'FontSize',fontsize)
xlim([0 error_max])

set(gcf,'PaperPositionMode','auto');
saveas(gcf, 'figure2_1b.png');

        % B3-Bipartisan Control;
figure('units','normalized','position',[.1 .1 0.5 0.2]);
subplot(1,2,1);
plot(grid, b_pcpossig(:,3),'k', 'LineWidth',width_line);
ylabel('Proportion Positive and Significant', 'FontSize',fontsize)
ylim([0 ylim_up])
xlabel('Standard Deviation of Error', 'FontSize',fontsize)
xlim([0 error_max])

subplot(1,2,2);
plot(grid, b_pcnegsig(:,3), 'k', 'LineWidth',width_line);
ylabel('Proportion Negative and Significant', 'FontSize',fontsize)
ylim([0 ylim_up])
xlabel('Standard Deviation of Error', 'FontSize',fontsize)
xlim([0 error_max])

set(gcf,'PaperPositionMode','auto');
saveas(gcf, 'figure2_1c.png');

        % B4-Fraction- Unified Republican Control Coefficient;
figure('units','normalized', 'position',[.1 .1 0.5 0.2]);
subplot(1,2,1);
plot(grid, b_pcpossig(:,4),'k', 'LineWidth',width_line);
ylabel('Proportion Positive and Significant', 'FontSize',fontsize)
ylim([0 ylim_up])
xlabel('Standard Deviation of Error', 'FontSize',fontsize)
xlim([0 error_max])

subplot(1,2,2);
plot(grid, b_pcnegsig(:,4),'k', 'LineWidth',width_line);
ylabel('Proportion Negative and Significant', 'FontSize',fontsize)
ylim([0 ylim_up])
xlabel('Standard Deviation of Error', 'FontSize',fontsize)
xlim([0 error_max])

set(gcf,'PaperPositionMode','auto');
saveas(gcf, 'figure2_1d.png');

        % B5-Fraction-  Bipartisan Control;
figure('units','normalized','position',[.1 .1 0.5 0.2]);
subplot(1,2,1);
plot(grid, b_pcpossig(:,5),'k', 'LineWidth',width_line);
ylabel('Proportion Positive and Significant', 'FontSize',fontsize)
ylim([0 ylim_up])
xlabel('Standard Deviation of Error', 'FontSize',fontsize)
xlim([0 error_max])

subplot(1,2,2);
plot(grid, b_pcnegsig(:,5),'k', 'LineWidth',width_line);
ylabel('Proportion Negative and Significant', 'FontSize',fontsize)
ylim([0 ylim_up])
xlabel('Standard Deviation of Error', 'FontSize',fontsize)
xlim([0 error_max])

set(gcf,'PaperPositionMode','auto');
saveas(gcf, 'figure2_1e.png');




%% Part 2.2. Do Monte Carlo simulation for SM (2015) with DW SE- procedure as before;

k=size(X_sm,2);
b_smse=zeros(k,rep);
se_smse=zeros(k,rep);
tstat_smse=zeros(k,rep);


rng(1991);
 for i=1:rep;
     dw_extended_noise=dw_extended_112+normrnd(0,dw_extended_se);
     
     [libs_senoise,libslag_senoise]=create_liberal(sm_congress,sm_state,dw_congress_ext,dw_state_ext,dw_extended_noise);
     
     Y_noise=libs_senoise;
     X_noise=[stfracmm, repgerr,bipartgerr, interac_rep,interac_bi, votlib2,blackpop,hispop,libslag_senoise,c];
   
     [b_smse(:,i), se_smse(:,i),tstat_smse(:,i)]=ols_weighted(X_noise,Y_noise,w_sm);
 end


    % Table with original SM extended model regression and the share
    % significant,pos and neg;
tstat_sig_mean=mean((abs(tstat_smse')>t_critical)*100)';
tstat_pos_mean=mean((tstat_smse'>t_critical)*100)';
tstat_neg_mean=mean((tstat_smse'<-t_critical)*100)';

    % Create significance stars for table display;
star_smse=cellstr(' ');

for i=1:k;
    if abs(tstat_sm(i))> norminv(0.95,0,1);
        star_smse(i)=cellstr('*');
    end
    
    if abs(tstat_sm(i))> norminv(0.975,0,1);
        star_smse(i)=cellstr('**');
    end
    
    if abs(tstat_sm(i))> norminv(0.995,0,1);
        star_smse(i)=cellstr('***');
    end
end
star_smse=star_smse';
    
    % Create table;
    
rownames_table_sm={'Fraction Majority-Minority'; 'Unified Republican Control'; 'Bipartisan Control'; 'Fraction x Republican'; 'Fraction x Bipartisan';
        'Voter Liberalism'; '\% Black'; '\% Latino'; 'Lagged Left of Median'; 'Intercept'};
SMCoeff=strcat(num2str(round(b_sm,4)),star_smse);
SMSE=round(se_sm,4);

table2_2=table(SMCoeff, SMSE , tstat_sig_mean, tstat_pos_mean, tstat_neg_mean, 'RowNames', rownames_table_sm);
writetable(table2_2,'table2_2.csv','Delimiter',',', 'WriteRowNames',true)

